function [new_node, new_elem] = process_grid_interpolation(node, elem, basis_type)
switch basis_type
    case "P1b"
        P = zeros(size(elem,1),2);
        for k = 1:size(elem,1)
            V = [node(elem(k,:),:)];
            P(k,:) = polygon_centroid(V);
        end
        new_node = [node; P];
        new_col = (Np+1:Np+size(elem,1))';
        new_elem = [elem, new_col];
    case "P2"
        % v0 = reshape(elem', 3*size(elem,1), 1);
        % v1 = reshape(([elem(:,2:3), elem(:,1)])', 3*size(elem,1), 1);
        % totalEdge = [v0,v1];
        % k = 1;
        % index1 = [];
        % for i = 1:size(totalEdge,1)
        %     for j = i+1:size(totalEdge,1)
        %         if all(totalEdge(i,:) == flip(totalEdge(j,:))) && all(totalEdge(i,:) ~= [0, 0])
        %             totalEdge(j,:) = [0, 0];
        %             index1(k,:) = [j,i];
        %             k = k+1;
        %         end
        %     end
        % end
        % uniqueEdge = totalEdge(~(totalEdge(:,1) == 0),:);
        % P = (node(uniqueEdge(:,1),:) + node(uniqueEdge(:,2),:))/2;
        % index1 = sortrows(index1);
        % index2 = [find(totalEdge(:,1)~=0), (1:size(P,1))'];
        % new_col = zeros(size(elem,1),3);
        % k = 1;
        % for i = 1:3:size(totalEdge,1)-2
        %     new_col(k,:) = [i, i+1, i+2];
        %     k = k+1;
        % end
        % for k = 1:size(index1,1)
        %     for i = 1:size(elem,1)
        %         for j = 1:3
        %             if new_col(i,j) == index1(k,1)
        %                 new_col(i,j) = index1(k,2);
        %             end
        %         end
        %     end
        % end
        % for k = 1:size(index2,1)
        %     for i = 1:size(elem,1)
        %         for j = 1:3
        %             if new_col(i,j) == index2(k,1)
        %                 new_col(i,j) = index2(k,2);
        %             end
        %         end
        %     end
        % end
        % new_node = [node; P];
        % new_elem = [elem, new_col+size(node,1)];
        v0 = reshape(elem', 3*size(elem,1), 1);
        v1 = reshape(([elem(:,2:3), elem(:,1)])', 3*size(elem,1), 1);
        totalEdge = [v0,v1];
        index1 = zeros(size(totalEdge,1),2);
        k = 1;
        for i = 1:size(totalEdge,1)
            if ~all(totalEdge(i,:) == [0, 0])
                for j = i+1:size(totalEdge,1)
                    if all(totalEdge(i,:) == flip(totalEdge(j,:)))
                        totalEdge(j,:) = [0, 0];
                        index1(k,:) = [j, i];
                        k = k+1;
                    end
                end
            end
        end
        index1 = index1(~(index1(:,1) == 0),:);
        uniqueEdge = totalEdge(~(totalEdge(:,1) == 0),:);
        P = (node(uniqueEdge(:,1),:) + node(uniqueEdge(:,2),:))/2;
        index1 = sortrows(index1);
        index2 = [find(totalEdge(:,1) ~= 0), (1:size(P,1))'];
        new_col = zeros(size(elem,1), 3);
        k = 1;
        for i = 1:3:size(totalEdge,1)-2
            new_col(k,:) = [i, i+1, i+2];
            k = k+1;
        end
        for k = 1:size(index1,1)
            idx = (new_col == index1(k,1));
            new_col(idx) = index1(k,2);
        end
        for k = 1:size(index2,1)
            idx = (new_col == index2(k,1));
            new_col(idx) = index2(k,2);
        end
        new_node = [node; P];
        new_elem = [elem, new_col+size(node,1)];
    otherwise
        error("Invalid basis type.");
end
end

%% polygon_centroid
function center = polygon_centroid(V)
n = size(V,1);
A = 0;
center_x = 0;
center_y = 0;
for i = 1:n
    j = mod(i,n) + 1;
    xi = V(i,1);
    yi = V(i,2);
    xj = V(j,1);
    yj = V(j,2);
    factor = xi * yj - xj * yi;
    A = A + factor;
    center_x = center_x + (xi + xj) * factor;
    center_y = center_y + (yi + yj) * factor;
end
center_x = center_x / (3 * A);
center_y = center_y / (3 * A);
center = [center_x, center_y];
end